The shrinking instability of toroidal liquid droplets in the Stokes flow regime 
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We analyze the stability and dynamics of toroidal liquid droplets. In addition to the Rayleigh 
instabilities akin to those of a cylindrical droplet there is a shrinking instability that is unique to 
the topology of the torus and dominates in the limit that the aspect ratio is near one (fat tori). We 
first find an analytic expression for the pressure distribution inside the droplet. We then determine 
the velocity field in the bulk fluid, in the Stokes flow regime, by solving the biharmonic equation for 
the stream function. The flow pattern in the external fluid is analyzed qualitatively by exploiting 
symmetries. This elucidates the detailed nature of the shrinking mode and the swelling of the cross- 
section following from incompressibility. Finally the shrinking rate of fat toroidal droplets is derived 
by energy conservation. 

PACS numbers: 47.20.Dr 
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Liquid droplets of various shapes are ubiquitous in na- 
ture pure and applied. They are found in rain, clouds, 
paint, lubricants, inks, dyes and oil [TIB and are being in- 
creasingly exploited in microfluidics [3] . The instabilities 
of liquid droplets have attracted attention since the be- 
ginning of the 19th century Early work of Plateau 
showed that a long cylindrical liquid droplet is unstable 
to capillary wave deformations of wavelength exceeding 
the droplet circumference. Rayleigh subsequently de- 
termined the most unstable capillary mode by solving 
the Navier- Stokes equation [TO, El- Purely planar liquid 
droplets are, in contrast, stable since capillary waves al- 
ways increase the droplet surface area and hence the free 
energy [l2|. Droplet instabilities thus probe the com- 
bined influence of the three-dimensional geometry of the 
droplets and their surface tension Q. 

In this paper we study the instability of liquid droplets 
in the form of three-dimensional axially symmetric solid 
tori, inspired by recent experiments in which bulk liquid 
tori are created by extruding water or glycerin through 
a metallic needle into a rotating bath of viscous silicone 
oil [l3| . Thin toroidal droplets exhibit Rayleigh instabil- 
ities analogous to those of the cylinder [10|, uA, UM •> w ^ n 
the additional requirement that the most unstable mode 
has wavelength A c commensurate with the outer circum- 
ference of the torus. When the outer circumference is 
an integer (n) times A c , the toroidal droplet eventually 
fissions into n solid spherical droplets (three-dimensional 
balls). Thus the change in topology of the droplet (solid 
torus breaks up into n balls) is governed by a Bohr quan- 
tization condition with the final number of balls playing 
the role of the principal quantum number n. Toroidal 
droplets also exhibit a fundamentally different type of 
instability in which the torus shrinks to close its interior 
hole, eventually becoming a single ball. This instability 
is a signature of the topological character of the torus 
and does not exist for a cylinder. Although it is present 
for a torus of any aspect ratio, it is preempted by the 
Rayleigh instability unless the torus is sufficiently fat (see 
Appendix B). 

The outline of this paper is as follows: we first analyze 
the shrinking instability of toroidal droplets by minimiz- 



ing a free energy controlled by interfacial surface ten- 
sion. We then derive the pressure distribution driving 
bulk flow of the fluid. The shrinking mode is then exam- 
ined in more depth via the Stokes equation, which is the 
large Ohnesorge number (Oh = rj/y/ pLcr » 1) limit of 
the Navier- Stokes equation. The biharmonic equation for 
the stream function determines the velocity field inside 
the toroidal liquid droplet. The shrinking of the droplet 
and simultaneous swelling of the cross-section from vol- 
ume conservation are clearly revealed in the flow. Finally, 
we calculate the shrinking speed by balancing the rate of 
free energy gain with the viscous dissipation rate. 




FIG. 1: Schematic of coordinates for a solid torus generated 
by rotating a circular disk of radius i?2, centered a distance 
Ri from the origin, around the vertical (z) axis. 

A three-dimensional axially symmetric solid torus is 
characterized by coordinates {u 1 = a,u 2 = 6,u 3 = r}, 
as shown in FigJTJ Here a is the angle around the tube, 
is the azimuthal angle around the z-axis, and r is the 
radial coordinate of the tube. The central circle of the 
solid torus with radius R\ (at r = 0) will be called the 
reference circle of the solid torus. The outer radius of 
the tube is denoted i?2- The aspect ratio of the toroidal 
surface is then <fi = R1/R2. The non-zero components of 
the metric tensor of the solid torus are gn = r 2 , #22 = 
(Ri -f rcosa) 2 and #33 = 1. 

The instability of a toroidal liquid droplet to shrinking 
can be seen in terms of the free energy F = crA. In 



2 



the shrinking process, R\ decreases and R 2 increases due 
to volume conservation. The change of free energy with 
radius is 



dF 
~dRi 



= 27r 2 aR 2 > 0. 



(1) 



Thus toroidal liquid droplets shrink to reduce the free en- 
ergy. Although this static analysis reveals the shrinking 
mode, the free energy alone does not provide a complete 
description of the system. In particular, determining the 
shrinking rate requires a study of droplet hydrodynamics. 

We first analyze the pressure distribution in a toroidal 
droplet to understand the driving force for bulk flow. 
Taking the divergence of the Navier- Stokes equation for 
an incompressible fluid shows that the pressure must be 
harmonic 



Ap(r, a) = 



(2) 



The boundary condition is given by the distribution of 
Laplace pressure on the interface between the inner and 
outer fluids 



p-po = crH, 
where H is the mean curvature H 



-R1+2-R2 cos a 



(3) 
For 



R2 (-R1+-R2 cos a) ' 

simplicity, we first consider the external pressure po as 
constant. The problem of solving for the pressure dis- 
tribution in the bulk fluid is then reduced to solving 
Laplace's equation, Eq.©, with the specified boundary 
condition. The Laplace pressure drop from the exterior 
(a = 0) to the interior (a = ir) of the torus is given by 
P(a = 0) — P(a = 7r) = 2<7 and is a measure of 
the asymmetry of the torus. Since (j) > 1, the Laplace 
pressure on the exterior of the toroid is always bigger 
than on the interior. One also sees that the asymmetry 
is more pronounced for a fat torus with aspect ratio (j) ap- 
proaching one. In the limit <fi —> 00, a toroid approaches 
a solid cylinder and the asymmetry as well as the shrink- 
ing mode disappear. Note that for the opposite case of 
constant pressure in the inner fluid, the pressure in the 
outer fluid will fall from the interior to the exterior of the 
torus. The outer fluid will therefore flow outward and 
the inner fluid will correspondingly flow inward, shrink- 
ing the droplet. Shrinking is thus a universal feature of 
one toroidal fluid inside another. 

Laplace's equation for the pressure Eq.(|2]) separates in 
the coordinates {p, ip, #}[l6|, [TtJ defined by 



a sinh p cos 6 
cosh p— cos ip 
a sinh p sin 9 
cosh p— cos (p 

a sin <£ 
cosh p— cos (p 



where a = r sinh p. 

Exploiting azimuthal symmetry and imposing the re- 
quirement that the pressure be finite as one approaches 
the reference circle, the physically acceptable solution 
takes the form [ItJ 
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a/ cosh p — cos Tp ^ a p Re[Q p _ 
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FIG. 2: The pressure p(r, a) vs. angle a at different radial 
distances away from the reference circle for R\ — 5 mm and 
aspect ratio <j> = R1/R2 = 1.5. The solid black curve is the 
pressure distribution on the boundary p = aH. Green curve: 
r — Ri/1.5 (boundary). Blue curve: r = Ri/3. Red curve: 
r = Ri/10. Dashed Black curve: r = i?i/100. The green 
curve fits the exact pressure on the boundary very well. 



where Re[Q p _i/ 2 (%)] is the real part of the associated 
Legendre function of the second kind. The coefficients 
a p can be determined by imposing the boundary con- 
dition in Eq.©. The pressure distribution inside the 
droplet is plotted in Fig|2] for aspect ratio (j) = 1.5, 
surface tension a = 4 x 10 -3 N/m and B\ = 5 mm. 
The pressure clearly drops from the exterior (a = 0, 2tt) 
to the interior (a = n). This pressure gradient drives 
the fluid towards the center of the toroid. As r de- 
creases the pressure distribution becomes more isotropic 
(^-independent). Note that only the p = mode in 
Eq.(|2j) contributes to the pressure near the reference cir- 
cle. This enables us to study the behavior of the fluid 
near the reference circle analytically. By inserting the 
zero mode in Eq.(j4j) into the Stokes equation one can 
determine the velocity field near the reference circle in 
z} coordinates: v x = v a sina — v r cosa = c, and 
v z — v a cosa + v r sina = 0. The velocity field near 
the reference circle is uniform towards the center of the 
toroid. 

We now turn to the velocity distribution in a viscous 
toroidal liquid droplet. In this regime of large Ohnesorge 
number the Navier-Stokes equation reduces to the Stokes 
equation 



Av = -Vp, 
V 



(5) 



where r\ is the fluid viscosity. Here viscous dissipation 
dominates over kinetic energy damping: 



d[\ J pv 2 dV]/dt pLcr 

OC ~z <C <C 1, 



/ dV<r' ij V i vi 



(6) 



where g[- 



r]{S/iVj + VjVi) is the viscous stress ten- 
sor. Since the characteristic speed of the fluid is much 
smaller than the speed of sound, we can treat the fluid as 
incompressible (V • v = 0)[18]. For incompressible fluids 
one can write the velocity field as the curl of a vector 
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FIG. 3: The velocity field inside and outside a cross-section 
of a toroidal liquid droplet. The dashed semi-circle is the 
interface of two distinct fluids. The velocity field inside the 
droplet is calculated by solving the biharmonic equation. The 
external flow pattern is schematically plotted by imposing 
boundary conditions and exploiting symmetry. Parameters 
are mode number v — 1, R\ — 5 mm and R2 = 2 mm. 



potential ip (v = V x ip) leading directly to 

aV= 0. 



(?) 



The complete velocity field can be obtained by solving 
the biharmonic vectorial equation Eq.© which reduces 
to a simplified scalar differential equation in the {p, (/?, 0} 
coordinates 191 



E 2 (E 2 ^j) = 0, 



(8) 



where ip, the stream function, is the only non-zero com- 
ponent ipo of the vector potential ip and the second- 
order partial differential operator E is given by E 2 = 
wh 2 [d p {^dp) +#<p(^0<p)], with w = a sinh(p)/[cosh(p) - 
cos((/?)] and h = [cosh(p) — cos(y?)]/a. Imposing the phys- 
ical requirements that approaching the reference circle v x 
tends to a finite value and v z — >• (reflection symmetry) 
yields the complete solution 



a sinh p 



+00 



(cosh p - 



c u sm(is(p)Ql_ s/2 (coshp)(9) 



Note that v x — » — and v z —> as p —> 00 (r — )> 0). 
Thus only the u = 1 mode contributes to the flow near 
the reference circle. The coefficients c v in Eq.(j9]) can 
be determined by matching to the velocity field on the 
interface. Assuming that high viscosity fixes the fluid 
particles on the interface to move with the interface as it 
shrinks, the boundary conditions are found to be v x o = 
V(l — 2^7) an d v z0 — — where xq and zq denote 
spatial points on the boundary and V = dR\/dt. The 
point (x = 0, z = 0) is the center of the cross section. The 
v = 1 mode is sufficient to fit this boundary condition. 

The velocity field inside the droplet in the laboratory 
frame is plotted in Figj3j The shrinking of the droplet 
is clearly indicated by the inward directed flow inside 
the droplet. One also sees that outer fluid within the 
toroidal hole is squeezed out. Further insight is provided 
by plotting the velocity field (see Fig|4]) inside the droplet 
in a reference frame comoving with the shrinking of the 
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FIG. 4: The velocity field inside the toroidal droplet of Fig[3] 
in a comoving reference frame shrinking with the droplet. The 
swelling of the cross section is readily inferred. 




FIG. 5: (a) The vorticity Q,e(r)/Q.e(r — R2) versus r. Param- 
eters: v = 1, ci = 1, a = 1 rad, Ri = 5 mm and R2 = 4mm . 
The vorticity field falls to zero at the reference circle, (b) Qo 
versus angle a for modes v in the range (—3,3). The vortic- 
ity Qo is rescaled to show different modes in the same figure. 
Parameters: (f) = 5, a = 1. Blue curve: v — (vorticity van- 
ishes). Black curve: v = 1. Brown curve: v = 2. Gray curve: 
v — 3. Green curve: v = — 1. Red curve: v — —2. Pink 
curve: v = — 3. Note that the number of peaks or valleys is 
determined by the mode number v. 



droplet. The swelling of the cross section resulting from 
volume conservation is clearly visible. 

The vorticity field, Q = V x v = —Aip is plotted in 
Fig|5] which shows its only non-zero component Q$ as 
a function of r and a respectively. FigJSfa) shows that 
the vorticity field is only significant near the boundary - 
it decays rapidly as one approaches the reference circle. 
FigEJb) shows Qq versus a for v G [—3,3]. The vorticity 
field vanishes at a = and tt due to its odd parity. The 
sign of Qq reflects the chiral property of vortices. The 
number of peaks and valleys on the z > plane (i.e., 



a G [0, 7r]) is 



n = 



-i/, < 
1/ - 1, i/ > 1 



1, 



and is therefore completely determined by the mode v. 

In the process of shrinking the free energy gained is 
dissipated in viscous damping. By equating the rate of 
change of the free energy, Eq.(pQ), to the viscous dissipa- 
tion rate we can obtain the shrinking speed. We focus 
here on the experimentally explored case of a low vis- 
cosity (r\i) inner fluid immersed in a viscous (rj ) outer 
bath [13j. In this case the dissipation occurs almost en- 
tirely in the outer fluid. Applying Stokes' equation for an 
incompressible fluid the dissipation rate can be separated 
into two parts 



E v 



dVcJijdiVj 



I 



dfiCr'ijVj 



dVvjdicr'ij 



(10) 



Here dfi is the i component of the area element of the in- 
terface. The first term is the heat flux on the fluid bound- 
ary and the second term is the dissipation rate inside the 
bulk fluid. The second term can be related to the vortic- 
ity: / Vjdid'^dV = rj J VjAvjdV = -rj J v • (V x Q)dV. 
Since the Reynolds number of the external fluid is very 
small in the experimental setup (Re « 10 -4 [13]), we 
may take the external flow to be as irrotational (vanish- 
ing vorticity) by recalling the experiment of flow through 
a cylindrical solid: an irrotational-rotational flow tran- 
sition occurs at Re ~ 1, below which the flow is irrota- 
tional [20|| . In the shrinking process, the toroidal droplet 
moves in the external fluid. This is equivalent to flow 
through the toroidal droplet. Since the viscosity of the 
internal fluid is very small in this case, the internal dissi- 
pation can be neglected. Thus we need to calculate only 
the surface integral in Eq.(10) to obtain the dissipation 
rate. We need the viscous stress on the boundary to 
evaluate the surface integral. Rotational symmetry and 
the limiting condition rji/r) <C 1 imply both a' rQ and (j' ra 
vanish at the interface. To determine o' rr = 2r] d r v r we 
need the gradient of v r at the interface. Assuming that 
the fluid particles near the boundary move together with 
the interface during shrinking, as a result of the viscous 
external fluid, we have 



„ 1 _ ggi 5 

O r V r \i n f er f ace — 2 -tt2COSCy. 

-Tto 



(11) 



Inserting Eq.(ll) along with the velocity field on the in- 
terface into the surface term of Eq.(10) yields 

E vis = -24A ((^) 2 - l)R2RiR2 • (12) 

By equating the rate of change of the free energy from 
Eq.(pQ) and the dissipation rate Eq. ([T2|) , we have 

1 



12 (j) 2 (t) - 1/2 



(13) 




FIG. 6: The time evolution of the normalized inner droplet 
radius Ri n (t) / Rin(0) for initial aspect ratio 1.4. The the- 
oretical result is the solid curve and the experimental data 
points are taken from Ref.[13l]. Parameters: Ri(0)= 3 mm; 
v s — cr/rjo — 133 /xm/s. 



and the interior hole of the droplet decreases in size ac- 
cording to 



Rin{t) — 



v s 2 <j>(t) + 1 
12 2 (t) - 1/2' 



(14) 



where (j)(t) = R\ (t) / R2 (t) and v s = a/rj Q . The shrinking 
speed is controlled by the aspect ratio of the droplet, 
in accord with experimental observations [T^j. In the 
limit of infinite aspect ratio (the cylinder) the shrinking 
speed vanishes, as required. The constant 1/2 in the 
denominator of Eq. (fT3|) plays an important role in the 
limit that the aspect ratio approaches one (fat tori) . The 
plot of Ri n (t) / Ri n (0) versus t is shown in Fig J6] for an 
initial aspect ratio R 1 (0)/R 2 (0) = 1.4, rji/rj = 1/30,000, 
<j = 4 mN/m, R±(0) ~ 3 mm and r} % = 10~ 3 kg/ (m • 
s). For these parameters v s = cr/r] w 133 /xm/s. Figj6] 
shows that droplets shrink with roughly constant speed, 
as found in [13[. Our results predict that it would take 
about 21 s (aspect ratio 1.4) and 50 s (aspect ratio 1.9) 
for a toroidal droplet to shrink to close the inner hole 
of toroid, in qualitative agreement with the experimental 
values of 25 s and 38 s respectively. Thus thinner toroidal 
droplets shrink more slowly, consistent with experimental 
observations [l3| . 

Our energy conservation approach to determining the 
shrinking speed can also be applied to a 2-dimensional 
system where it yields an analytical result. Consider 
a shrinking hole on a liquid film. The limiting case of 
a shrinking toroidal liquid droplet with R in — ^ and 
Vi/^o » 1 can be modelled as such a 2-dimensional 
system, since the dynamics of the hole becomes indepen- 
dent of the fluid far away from the hole. As the hole 
shrinks, a flow will be induced outside the hole on the 
film. In the Stokes flow regime, the velocity field can 
be derived analytically in the polar coordinates {p, 0} as 
v p = ^r(t)r(t) and vq = 0. By energy conservation, the 
shrinking speed of the hole can also be derived analyt- 
ically. By equating the rate of change of the line en- 
ergy E s = = 27T7 r and the viscous dissipation rate 
E V is = — J pdpd6(a fpp d p v p ) = —2i\r] r 2 , we have r = ^, 
where rj is the viscosity of fluid and 7 is the line tension. 
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We expect that the formalism employed here will have 
a variety of applications to the dynamics of fluid in- 
terfaces. It may also be extended to liquid crystalline 
droplets where the interplay of liquid crystalline order 
and the shape of the droplet should be very rich. 



For tori with typical aspect ratios R1/R2 = 10 and 2, 
we find the ground states in the {c2,C3|c2,C3 G [— 6, b]} 
space shown in FigJTl The cross-sections are very close 
to circular. In the experimental work of [l3| this is also 
true. 
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Appendix 

A. The shape of cross section 

So far we have assumed that droplets remain perfectly 
circular in cross-section as they shrink. Here we show 
that this assumption is well justified. 
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FIG. 7: Minimum energy cross-sectional shapes within a two 
parameter family of possible shapes (see Ea. (|A.l[) ). The 
dashed curves are the unperturbed circular shapes: (a) as- 
pect ratio 10 (thin torus) and (b) aspect ratio 2 (fat torus). 
Parameters: R2 = 1. 

The shape of a toroidal liquid droplet is characterized 
by the radii R\ and R2 which may in general vary with 
a and 0. Retaining azimuthal symmetry we consider the 
following ansatz for R2 at a fixed time: 



R 2 (a) 



C2P2(cosa) -f c3P3(cosa). 



(A.l) 



The second term describes an ellipse which is symmetric 
about z axis, while the third term describes a shape with 
three round corners, which is asymmetric about the z- 
axis (we are ignoring the shrinking mode here, described 
by a Pi (cos a) term). The shape of the droplet is specified 
by points in the {02,03^2,03 G [—&,&]} space. 

We numerically search for the ground state in the 
{c2,C3|c2,C3 G [—&,&]} space for which 



L = A- A + \\V-V 



(A.2) 



is minimized, where Vo and Aq are the volume and surface 
area of the unperturbed droplet. A is set to be large to 
impose volume conservation. We take a = 1 and b = 0.2. 



It is observed experimentally [l3j that the Rayleigh in- 
stability disappears for sufficiently fat solid tori (Ri (t = 
0)/i?2(£ = 0) < 2) whereas the shrinking mode is present 
for all aspect ratios. Here we derive a lower bound on the 
aspect ratio for the emergence of the Rayleigh instability. 

Two conditions must be satisfied for the Rayleigh in- 
stability: (1) modes with wavelength A > A c , where A c 
is the minimum wavelength of the Rayleigh instability 
mode and (2): 



u k (t = t 1 )>R 2 (t = 0), 



(B.l) 



where Uk is the perturbation amplitude and t\ is the life- 
time of the shrinking droplet (Ri n (t = t\) = 0). It is well 
known [TO, [H| that Uk grows exponentially: 



u k {t)=u k {<d)e Vet/R2( °\ 



(B.2) 



where v s = a/77 is the characteristic speed and i?2(0) the 
characteristic length scale of the system. We assume that 
Uk grows exponentially all the way until breakup of the 
droplet. On the other hand we have shown that Ri n (t) 
decreases almost linearly in time. So formally, we have 



R in (t) = R in (t = 0) - v c t, 



(B.3) 



from which we have t\ = Rtn ^ °) . By inserting 
Eq. (jB.2IB.3|) into Eq.flEU, we obtain 



fli(O) 



> 1 



cm — 7-f , 



(B.4) 



where c — v c /v s is an aspect ratio factor of order one 
that tends to 1/2 for aspect ratio one by Eq. (fT4|) . Thus 
the Rayleigh instability is dominant for sufficiently thin 
tori. 

It can be checked that for aspect ratios satisfying 
Eq. (|B.4[) even the perimeter of the interior of the torus 
(27r(i?i(0) — i?2(0))) can accommodate the Rayleigh in- 
stability mode, i.e., 



ifc(0) 



wfe(0) 
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